{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "import joblib\n",
    "from scipy.fftpack import fft\n",
    "from scipy.signal import find_peaks\n",
    "from scipy.stats import skew, kurtosis\n",
    "\n",
    "# 加载训练好的随机森林模型\n",
    "model_path = 'rf_model.pkl'  # 确保路径正确\n",
    "rf_model = joblib.load(model_path)\n",
    "\n",
    "file_path_train = r'C:\\Users\\Administrator\\Desktop\\聚合后的材料数据.xlsx'  # 修改为实际文件路径\n",
    "\n",
    "# 定义变量范围\n",
    "temperatures = [25, 50, 70, 90]\n",
    "materials = ['材料1', '材料2', '材料3', '材料4']\n",
    "waveforms = ['三角波', '梯形波', '正弦波']\n",
    "freq_min, freq_max = 49990, 501180\n",
    "\n",
    "# 读取 Excel 文件，从 F 列的第二行开始提取磁通密度数据\n",
    "flux_data = pd.read_excel(file_path_train).iloc[1:, 5:].values  # 从第二行、F列开始提取数据\n",
    "\n",
    "# 使用 np.nan_to_num 处理 NaN 和无穷大值，避免后续报错\n",
    "flux_data = np.nan_to_num(flux_data, nan=0.0, posinf=1e6, neginf=-1e6)\n",
    "\n",
    "# 1. 提取特征，包含峰值、时域和频域特征\n",
    "def extract_features(temp, freq, material, waveform, flux_data_row):\n",
    "    # 清理 flux_data_row 中的 NaN 和无穷大\n",
    "    flux_data_row = np.nan_to_num(flux_data_row, nan=0.0, posinf=1e6, neginf=-1e6)\n",
    "\n",
    "    # 检查 flux_data_row 是否全为零，防止异常\n",
    "    if np.all(flux_data_row == 0):\n",
    "        print(\"Invalid flux_data_row: All values are zero\")\n",
    "        return np.zeros(15)  # 返回默认特征，防止后续计算出错\n",
    "\n",
    "    # 提取磁通密度中的峰值（每行的最大值）\n",
    "    max_flux_density = np.max(flux_data_row)\n",
    "\n",
    "    # 提取时域特征\n",
    "    time_features = [\n",
    "        np.mean([temp, freq, max_flux_density]),  # 均值\n",
    "        np.std([temp, freq, max_flux_density]),  # 标准差\n",
    "        max([temp, freq, max_flux_density]),  # 最大值\n",
    "        min([temp, freq, max_flux_density]),  # 最小值\n",
    "        (max_flux_density - temp) / len([temp, freq, max_flux_density]),  # 斜率\n",
    "        len(find_peaks(flux_data_row)[0]),  # 波峰数\n",
    "        np.sum(np.abs(np.diff([temp, freq, max_flux_density]))),  # 平滑度\n",
    "        skew([temp, freq, max_flux_density]),  # 偏度\n",
    "        kurtosis([temp, freq, max_flux_density])  # 峰度\n",
    "    ]\n",
    "\n",
    "    # 提取频域特征\n",
    "    freq_features = [\n",
    "        np.argmax(np.abs(fft([temp, freq, max_flux_density]))),  # 主频率\n",
    "        np.sum(np.abs(fft([temp, freq, max_flux_density]))[1:10]) / np.sum(np.abs(fft([temp, freq, max_flux_density]))),\n",
    "        # 高频能量比\n",
    "        -np.sum((np.abs(fft([temp, freq, max_flux_density])) / np.sum(np.abs(fft([temp, freq, max_flux_density])))) *\n",
    "                np.log(np.abs(fft([temp, freq, max_flux_density])) / np.sum(np.abs(fft([temp, freq, max_flux_density]))))),\n",
    "        # 谱熵\n",
    "        np.sum(np.square(np.abs(fft([temp, freq, max_flux_density]))))  # 功率谱密度\n",
    "    ]\n",
    "\n",
    "    # 添加材料和波形的特征\n",
    "    material_feature = [materials.index(material)]  # 材料种类\n",
    "    waveform_feature = [waveforms.index(waveform)]  # 波形种类\n",
    "\n",
    "    # 合并时域、频域特征以及材料、波形特征，并包含从 flux_data_row 中提取的峰值\n",
    "    features = time_features + freq_features + material_feature + waveform_feature + [max_flux_density]\n",
    "\n",
    "    # 处理 NaN 或无穷值，防止模型预测错误\n",
    "    features = np.nan_to_num(features, nan=0.0, posinf=1e6, neginf=-1e6)\n",
    "\n",
    "    # 如果 features 长度不匹配，返回默认值以避免模型报错\n",
    "    if len(features) != 15:\n",
    "        print(f\"Invalid feature length: {len(features)}, expected 15.\")\n",
    "        return np.zeros(15)  # 返回默认值\n",
    "\n",
    "    return np.array(features).reshape(1, -1)\n",
    "\n",
    "# 2. 目标函数\n",
    "def objective_function(temp_idx, mat_idx, wave_idx, freq, flux_data_row):\n",
    "    try:\n",
    "        # 确保索引在有效范围内\n",
    "        temp_idx = max(0, min(int(round(temp_idx)), len(temperatures) - 1))\n",
    "        mat_idx = max(0, min(int(round(mat_idx)), len(materials) - 1))\n",
    "        wave_idx = max(0, min(int(round(wave_idx)), len(waveforms) - 1))\n",
    "\n",
    "        # 提取特征并进行预测\n",
    "        input_features = extract_features(temperatures[temp_idx], freq, materials[mat_idx], waveforms[wave_idx], flux_data_row)\n",
    "        input_features = input_features.reshape(-1)  # 确保输入为一维数组\n",
    "\n",
    "        if np.all(input_features == 0):  # 防止异常输入进入模型\n",
    "            print(\"All-zero input features detected, skipping this iteration.\")\n",
    "            return float('inf')\n",
    "\n",
    "        core_loss = rf_model.predict(input_features.reshape(1, -1))[0]  # 预测磁芯损耗\n",
    "\n",
    "        # 传输磁能是频率和磁通密度峰值的乘积\n",
    "        transmitted_energy = freq * np.max(flux_data_row)\n",
    "\n",
    "        # 定义目标函数：最小化磁芯损耗，并最大化传输磁能\n",
    "        epsilon = 1e-6\n",
    "        objective = core_loss / (transmitted_energy + epsilon)\n",
    "\n",
    "        # 检查 objective 值的有效性\n",
    "        if np.isnan(objective) or np.isinf(objective):\n",
    "            print(f\"Invalid objective value: {objective}\")\n",
    "            return float('inf')\n",
    "\n",
    "        return objective\n",
    "\n",
    "    except Exception as e:\n",
    "        print(f\"Error in objective function: {e}\")\n",
    "        return float('inf')\n",
    "\n",
    "# 3. 初始化种群\n",
    "def initialize_population(pop_size):\n",
    "    population = []\n",
    "    for _ in range(pop_size):\n",
    "        temp_idx = random.randint(0, len(temperatures) - 1)  # 温度索引\n",
    "        mat_idx = random.randint(0, len(materials) - 1)  # 材料索引\n",
    "        wave_idx = random.randint(0, len(waveforms) - 1)  # 波形索引\n",
    "        freq = random.uniform(freq_min, freq_max)  # 频率\n",
    "        particle = [temp_idx, mat_idx, wave_idx, freq]\n",
    "        population.append(particle)\n",
    "    return population\n",
    "\n",
    "# 4. PSO 更新函数\n",
    "def update_particle(particle, pbest, gbest, velocity, w, c1, c2):\n",
    "    new_velocity = [\n",
    "        w * velocity[i] + c1 * random.random() * (pbest[i] - particle[i]) + c2 * random.random() * (gbest[i] - particle[i])\n",
    "        for i in range(len(particle))\n",
    "    ]\n",
    "    new_particle = [particle[i] + new_velocity[i] for i in range(len(particle))]\n",
    "    return new_particle, new_velocity\n",
    "\n",
    "# 5. 优化算法函数 (PSO-GWO-DE)\n",
    "def pso_gwo_de_optimization(pop_size, max_iter, w, c1, c2, F, CR, flux_data):\n",
    "    population = initialize_population(pop_size)\n",
    "    velocity = [[random.uniform(-1, 1) for _ in range(4)] for _ in range(pop_size)]  # 初始化速度\n",
    "\n",
    "    pbest = population.copy()  # 个体最优\n",
    "    gbest = min(population, key=lambda x: objective_function(*x, flux_data[np.random.randint(len(flux_data))]))  # 全局最优\n",
    "    best_value = objective_function(*gbest, flux_data[np.random.randint(len(flux_data))])\n",
    "\n",
    "    values = []  # 记录目标函数值\n",
    "\n",
    "    for iteration in range(max_iter):\n",
    "        for i in range(pop_size):\n",
    "            # PSO 更新\n",
    "            new_particle, new_velocity = update_particle(population[i], pbest[i], gbest, velocity[i], w, c1, c2)\n",
    "            velocity[i] = new_velocity\n",
    "            if objective_function(*new_particle, flux_data[np.random.randint(len(flux_data))]) < objective_function(*pbest[i], flux_data[np.random.randint(len(flux_data))]):\n",
    "                pbest[i] = new_particle\n",
    "                if objective_function(*pbest[i], flux_data[np.random.randint(len(flux_data))]) < objective_function(*gbest, flux_data[np.random.randint(len(flux_data))]):\n",
    "                    gbest = pbest[i]\n",
    "\n",
    "        # 记录最\n",
    "        # 记录最优值\n",
    "        best_value = objective_function(*gbest, flux_data[np.random.randint(len(flux_data))])\n",
    "        values.append(best_value)\n",
    "        print(f\"Iteration {iteration + 1}, Best Value: {best_value}\")\n",
    "\n",
    "    return gbest, best_value, values\n",
    "\n",
    "# 6. 运行优化算法\n",
    "pop_size = 50\n",
    "max_iter = 1000\n",
    "w = 0.7\n",
    "c1, c2 = 1.5, 1.5\n",
    "F = 0.8\n",
    "CR = 0.9\n",
    "\n",
    "best_solution, best_value, values = pso_gwo_de_optimization(pop_size, max_iter, w, c1, c2, F, CR, flux_data)\n",
    "\n",
    "# 输出最优解\n",
    "optimal_conditions = {\n",
    "    '温度，oC': temperatures[int(best_solution[0])],\n",
    "    '磁芯材料': materials[int(best_solution[1])],\n",
    "    '励磁波形': waveforms[int(best_solution[2])],\n",
    "    '频率，Hz': best_solution[3],\n",
    "}\n",
    "\n",
    "print(\"最优条件：\")\n",
    "for key, value in optimal_conditions.items():\n",
    "    print(f\"{key}: {value}\")\n",
    "\n",
    "print(f\"最优目标函数值: {best_value}\")\n",
    "\n",
    "# 绘制优化过程的收敛图\n",
    "plt.figure(figsize=(10, 6))\n",
    "plt.plot(values, marker='o')\n",
    "plt.title('PSO-GWO-DE混合自适应算法优化过程')\n",
    "plt.xlabel('迭代次数')\n",
    "plt.ylabel('目标函数值')\n",
    "plt.grid(True)\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
